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Abstract 

It has been known for some time that human autoimmune diseases can be triggered by viral 
infections. Several possible mechanisms of interactions between a virus and immune system 
have been analysed, with a prevailing opinion being that the onset of autoimmunity can in many 
cases be attributed to "molecular mimicry", where linear peptide epitopes, processed from viral 
proteins, mimic normal host self proteins, thus leading to a cross-reaction of immune response 
against virus with host cells. In this paper we present a mathematical model for the dynamics 
of an immune response to a viral infection and autoimmunity, which takes into account T cells 
with different activation thresholds. We show how the infection can be cleared by the immune 
system, as well as how it can lead to a chronic infection or recurrent infection with relapses and 
remissions. Numerical simulations of the model are performed to illustrate various dynamical 
regimes, as well as to analyse the potential impact of treatment of autoimmune disease in the 
chronic and recurrent states. The results provide good qualitative agreement with available data 
on immune responses to viral infections and progression of autoimmune diseases. 

1 Introduction 

A successful immune system relies on its ability to discriminate between cells infected with a 
pathogen such as virus and the uninfected cells of the host. The breakdown in discrimination of 
self-antigens results in autoimmunity, where the immune system attacks specific cells or organs. 
The specificity of the immune response focuses disease within different target organs, for example, 
pancreatic /3-cells in insulin-dependent diabetes mellitus type-1 (IDDM), the central nervous system 
in multiple sclerosis (MS), or the retina in uveitis - an inflammation of the eye |44 [ 142 } 125]. 

Multiple factors are known to contribute toward the onset and development of autoimmune 
diseases, including genetic predisposition, age, and environment. Amongst the environmental fac- 
tors, the major identified triggers of autoimmunity are believed to be infectious pathogens [6j[9]. 
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Experimentally, direct infection of islet cells in the pancreas led to bystander damage of the islet 
cells and autoimmunity, and it was concluded that the release of sequestered antigen was critical to 
the development of disease [22]. There is also a strong association between infection with hepatitis 
C and autoimmune hepatitis [32J. While it is certainly not the case that all autoimmune disease is 
triggered by infection of the organ that subsequently becomes the target, this is a reasonable model 
for several known examples of organ specific autoimmunity. 

There are many mechanisms by which host infection by a pathogen may contribute to autoimmu- 
nity, including triggering the innate immune system, molecular mimicry and bystander activation 
(see, e.g., [TSj for a recent review). Molecular mimicry is thought to be particularly important 
when viruses cause autoimmunity [HI [T71 [331 EI]- An immune response elicited against a viral 
protein that mimics a self antigen will not only eliminate the virus, but can also target normal host 
cells that display the cross-reactive self-antigen. Such interactions may play a role both in causing 
disease and also in precipitating a relapse of disease. Between such triggering events, disease can 
of remit in part because of the existence of regulatory pathways that dampen immune responses. 

Whilst significant advances have been made in mathematical modelling of various aspects of 
general virus dynamics and the interactions between viruses and the immune system of the host, 
theoretical studies of autoimmunity have been quite limited in comparison. One of the earliest 
mathematical models of autoimmunity did not explicitly include specific causes of autoimmunity 
but rather concentrated on the interactions between effector and regulator cells, and used this to 
get an insight into T cell vaccination [36]. Further progress in the context of T cell vaccination 
was made by Borghans et al. [3 [8] , who demonstrated how the interactions of autoreactive and 
regulatory T cells can lead to the onset of autoimmunity or stable oscillations around a vaccinated 
state. Burroughs et al. [10] have studied the onset of autoimmunity through bystander activation. 
Wodarz and Jansen [52J analysed autoimmunity in the context of viral causes of cancer. They 
included viral infections indirectly through an increased rate of uptake of self-antigen by antigen- 
presenting cells. Leon et al. [28j EH [30] have studied the dynamics of interactions between different 
T cells for the purposes of regulation of immune response and control of autoimmune reaction. More 
recently, Iwami et al. |23[ [24] derived and studied a model for autoimmunity, which makes explicit 
account of the virus dynamics and its interaction with the immune system by means of linear or 
nonlinear immune response. Despite its simplicity, this model appears unable to reproduce a normal 
clearance of virus during a single infection, as it does not allow for a viral expansion. Various roles 
played by the regulatory T cells in the dynamics of autoimmunity have recently been analysed by 
Alexander & Wahl [2]. 

There are several ways to account for the ability of T cells to discriminate between cells pre- 
senting self antigens and infected cells. One of these is through regulatory cells, which are triggered 
by autoantigens and inhibit the activity of autoreactive T cells. This approach has already been 
extensively used in models of immune response, see Alexander & Wahl [2], Burroughs et al. |10j 
etc. Another approach is to consider T cells, which can perform a wide range of immune function 
by virtue of having different or tunable activation thresholds (TAT). This concept has been defined 
previously as: "activation is a threshold phenomenon and the threshold is tuned by the stimula- 
tory experience of the cell" |18j . i.e. T cells continually tune their responsiveness to T cell antigen 
receptor (TCR) stimulation through stimuli evoked by autoantigens. Because the degree of au- 
toreactivity of T cells is continuously controlled through their activation and tuning, this approach 
provides another way of modelling the onset and development of autoimmune disease. Grossman 
and Paul [18[ 120]. and Grossman and Singer |19j developed models with tunable activation thresh- 
olds that were applied to peripheral and to central T cell activation. Altan-Bonnet and Germain 
have modelled signalling threshold and shown differences in activation/response threshold that are 
dependent on the activation state of the T cell pQ. Noest [37] has shown how the need for activa- 
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tion threshold tuning arises from the first principles of signal detection theory, see also Scherer et 
al. |45j for further discussion of this issue, van den Berg and Rand have studied mathematically 
two cellular response models of the dynamics of tunable activation threshold |49| . Carneiro et al. 

have performed a comparative study of two mechanisms of self-tolerance: tuning of activation 
thresholds and control by specific regulatory T cells. The authors have shown that these two mech- 
anisms are complementary and together provide a plausible explanation of the observed dynamics 
of immune tolerance. Besides purely theoretical studies, dynamical changes in T cell activation 
during their circulation have also been shown experimentally both in the mouse and in man, where 
it has important implications for the outcome of specific therapeutic interventions [5] I3UI 14"3"1 14"5]. 

In this paper we propose and study a mathematical model for autoimmune disease caused 
by viral infections through molecular mimicry. By introducing separate populations of regular 
activated T-cells and T-cells with a lower activation threshold to self-antigen which arises as a result 
of the infection, the model is able to qualitatively reproduce normal aspects of immune behaviour. 
The model presented in this paper differs from an earlier work in this particular modelling aspect, 
as well as in the explicit account for a viral infection represented both by infected cells, and by a 
separate population of free virions. 

The organization of this paper is as follows. In the next section we discuss various biological 
assumptions behind the origins of autoimmunity and derive the corresponding mathematical model. 
In Section 3 we perform a systematic analysis of the steady states of the model and their stability. 
Section 4 contains results of numerical simulations in different parameter regimes, which illustrate 
clearance of infection, onset of autoimmunity, periodic flare-ups, as well as the dynamics of multiple 
infections. The paper concludes with the discussion of obtained results in Section 5. 

2 The model 

We are interested in modelling the interactions between a viral infection and a human host immune 
system, that lead to the onset of autoimmunity through molecular mimicry. When cells in a 
particular organ of the body become infected with a pathogen, the elicited immune response may 
cross-react with one or more self-antigens that share determinants with the pathogen [51]. This 
subsequently can lead to a breakdown of tolerance for self antigens through the appearance of 
lymphocytes capable of attacking the host's own cells, both in the same and in other organs. For 
both viral and bacterial pathogens, there is evidence that their antigens derived from pathogens 
engage and expand both CD4 and CD8 T cells that drive autoimmune disease [121 140j . Furthermore, 
recent genetic studies of autoimmunity implicate many genes that control T cell expansion and 
activation threshold [14J. Hence, it is important to understand a model of autoimmune disease 
that specifically addresses the dynamics of T cells. While it is unarguable that antibodies are 
important in many types of autoimmune disease, there is good evidence that T cell recognition 
of antigenic peptides is often a critical initiating step. For example, in a model of rheumatoid 
arthritis, antibodies were sufficient to induce disease [27] but the development of antibodies depends 
on prior T cell interactions with bacteria [55]. Furthermore, in some of the experimental models 
of autoimmunity, B cells are dispensable, and disease develops when they are not present [54]. 
Therefore, the balance of evidence is for a necessary T cell component in the onset and development 
of autoimmunity, and our model specifically seeks to address this. 

To model the dynamics of an immune response during a viral infection and possible onset 
of autoimmunity, we employ a model similar to those studied earlier in the context of immune 
responses [381 152j . Let A denote the number of susceptible cells in a particular organ or tissue. 
We consider two distinct populations of susceptible cells, A\ and Ai to allow for a situation when 
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autoimmunity takes place in a different organ to the one where the original infection occurs. It 
is assumed that in the absence of infection or autoimmunity, these two cell populations would be 
maintained at a certain constant level supported by homeostasis. It has been previously shown 
|23[ [23] that a specific form of the growth function for susceptible cells can have a significant effect 
on the overall dynamics of autoimmune diseases, and here we use the logistic form, as studied in 
Iwami et al. |23| . and by Perelson and Nelson |41] in models for HIV infection. 

When a person acquires a viral infection, a proportion of cells become infected with this virus, 
and after a certain period of time the infected cells, whose population is denoted by F, will start 
producing virions (free virus particles) V that will go on to infect other as yet uninfected cells A 2 
(and, possibly, A\). All time constants, such as the time to encounter an uninfected cell and a time 
required for cell entry, are implicitly included in the rate of infection A. In the absence of infection, 
the naive T-cells Tj n are taken to follow a logistic growth. Once activated, Tj n cells become T\ 
cells that have the ability to kill infected cells F. For a fraction of these cells T 2 their activation 
threshold for stimulation by susceptible cells that are not infected is reduced, allowing them to kill 
A cells. This is the autoimmune response. It has already been shown that T cells with varying 
activation threshold can have a significant effect on the dynamics of the immune response |18[ [19] . 

With the above assumptions, the model to be analysed in this paper has the form 

<,Al r 1 A 1 (l-4^ - P1XA1V - a a T 2 A h 



dt V Aq 



<Lh r 2 A 2 ( 1 - ^ ] - XA 2 V - P2 a a T 2 A 2 , 



dF_ 
~dt 



dt ' ' \ N 2/ 

\(p 1 A l + A 2 )V - (J-fF - apTxF - a a T 2 F, 
dT- / T 



g t T in ( 1 - ) - a act T in F, i 1 ) 



dt * \ M 

a act T in F - Ti(/ii + r), 

tTi - fi 2 T 2 , 



~dT 



dT 2 
dt 



dV 
~dt 



kF - 7 V, 



where A% !2 , F, Ti n , T\, T 2 and V denote the populations of susceptible cells (possibly, affected by 
infection and autoimmunity separately), infected cells, naive T cells, activated T cells, T cells with 
a lower activation threshold to self-antigen, and free virus, respectively. As it has already been 
mentioned, in the absence of infection which can trigger autoimmune reaction, the susceptible cells 
A± i2 reproduce logistically with their respective linear growth rates r\^ 2 and carrying capacities 
Aq^. These cells become infected at a rate A, and they are destroyed at a rate a a by autoreactive 
T-cells T 2 with a lower activation threshold to self-antigen (this is the actual implementation of the 
autoimmunity mechanism). Naive T-cells Ti n are assumed to have a logistic growth with a linear 
growth rate gt and a carrying capacity M; they get primed by dendritic cells at a rate taken to be 
proportional to the number of infected cells with a constant a ac t- 

Once infected, cells presenting foreign antigen F die at a rate fiF, and they are also destroyed 
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by the activated T-cells T\ at a rate ap, and by the autoreactive T-cells T2 at a rate a a . Activated 
T-cells T\ die at a rate fj,\, and at a rate r they produce autoreactive T cells T2 with a lower 
activation threshold to self-antigen. Finally, the T cells T 2 die at a rate [12- Free virions V are 
produced by an infected cell at a rate k, and 7 is the natural clearance rate of the virus. 

Parameters < pi,p% < 1 control whether the populations of susceptible cells are affected 
by autoimmune reaction and/or infection. When pi = P2 = 0, only cells A2 are affected by the 
infection, and the autoimmunity only affects the population Ai of other cells. Biologically, this 
situation arises when a pathogen causes infection in one part of the body, and the autoimmune 
response causes damage in other parts of the body. If p\ = and P2 > 0, only the cell population 
A2 experiences infection, but both of the susceptible cell populations are affected by the developing 
autoimmunity. Conversely, if P2 = and p\ > 0, both Ai and A2 cells become infected, and only 
Ai also experiences autoimmunity. Finally, in the case of p\ t 2 > 0, both of the cell populations 
A\ and A2 get exposed to both the infection and the autoimmune reaction, although they can 
potentially be affected by infection/autoimmunity at different rates, and the differences in intrinsic 
growth rates and carrying capacities of these two cell populations can also lead to a significantly 
different dynamics of these two cell populations. 

The main emphasis of the above model is on the separation of T cell populations into two 
activated populations, one of which is capable of an autoimmune reaction through having a lower 
activation threshold to self-antigen. As we are particularly interested in the role of foreign infections 
in the onset of autoimmunity, we have not explicitly included in the model several other aspects that 
can be of interest in specific contexts, such as antibody response or regulatory T cells. Furthermore, 
we have not included memory cells in our model, which can be important in the analysis of a longer- 
term dynamics or multiple infections. 

3 Steady states 

We begin our analysis of the system ([I]) by considering its steady states E* = (A*, A\, F* , T* n , T*, T£ , V*), 
which can be found by equating the right-hand side of system (TjQ) to zero. Stability is determined 
by the eigenvalues of the Jacobian of linearization of the system ([TJ near each of the steady states. 
In order to systematically study the steady states of system (1), we first consider all steady states 
with V* = (this immediately implies F* = T* = T% = 0). There are exactly eight of such steady 
states: 



E\ 


= (0,0,0,0,0,0,0), 


Et 


= (N 1: 


0,0,0,0,0,0), 


El 


= (0,0,0, A/, 0,0,0), 


Et 


= (N h 


0,0,M, 0,0,0), 


Et 


= (0,JV 2 , 0,0, 0,0,0), 


E$ 


= (N 1 . 


N 2 , 0,0, 0,0,0), 


E\ 


= (0,iV 2 ,0,M, 0,0,0), 


E* 


= (N h 


,N 2 ,0,M, 0,0,0) 



(2) 



Computing the Jacobian at each of these steady states shows that E^ to Ej are always saddles 
(and hence unstable) for any values of the system parameters. The steady state E$ is a stable 
node, provided 

k\(N 2 +piN 1 ) < pry, 

or a saddle if 

k\{N 2 +piN l ) > p, F j. 

At 

k\{N 2 +piN 1 ) = p Fl , 
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the steady state E% undergoes a steady state bifurcation, where one of the eigenvalues goes through 
zero along the real axis. Biologically, the above stability condition represents the rate of change in 
the number of infected cells being negative at the steady state. 

When V* ^ 0, we can have either A* = and A 2 / 0, or A\ / and A* 2 = 0, or A\ ^ and 
A2 / 0. In the case where A\ = and A 2 7^ 0, there are two options. One of these is Eg with 
Ti n = = T2 = 0, in which case we also have 

= M£7 F * = r _2l(i_ Vfl \ v* = — (l- Mf7 
2 Xk 1 Xk V XkN 2 J ' A V AfciV 2 

The steady state Eg with such values of the variables is stable, provided the following conditions 
are satisfied 

p 1 r 2 N 2 Xk > Xkr x N 2 + pir 2 jfi F , a> act r 2 jXkN 2 > g t X 2 k 2 N 2 + a act r 2 -/ 2 fi F , 

(3) 

7//Fr 2 (mf + 7) + XkN 2 [(/j, F + l) 2 + Mf7 - XkN 2 ] > 0. 

It is noteworthy that when p\ = 0, this steady state is unstable for any values of other parameters. 
For pi = 1, this steady state can only be stable, provided the linear growth rate r 2 sufficiently 
exceeds 77. The steady state Eg can undergo a steady state bifurcation when 

pir 2 N 2 Xk = XknN 2 +pir 2 ^ii F , or a act r 2 -fXkN 2 = g t X 2 k 2 N 2 + a act r 2 -f 2 fi F , (4) 

and it can also undergo a Hopf bifurcation, when the system parameters satisfy the conditions 

pir 2 N 2 Xk > XkriN 2 + pir 2 -ffi F , a act r 2 jXkN 2 > g t X 2 k 2 N 2 + a act r 2 j 2 fi F , 

(5) 

l^FriiliF + 7) + XkN 2 [(fj, F + 7) 2 + fi F ^y - XkN 2 ] = 0. 

An important observation is that the manifold A\ = Tj„ = T\ = T 2 = is flow- invariant, and 
the Hopf bifurcation of the steady state Eg takes place inside this manifold, which results in the 
appearance of a periodic orbit confined to the same manifold. The importance of this observation 
lies in the fact that this periodic solution does not cause oscillations of any of the variables which 
would go below zero, making it unrealistic. 

The second option for the case A\ = and A 2 / is given by the steady state E* , for which 
T* n satisfies the quadratic equation 



-r* \2 
-in) 



gt7Ct act [r 2 -f(a F fi 2 + a a r) + p 2 a a TN 2 Xk](T* 

-9t[r 2 l 2 a act M{a F ii 2 + a a r) - X 2 k 2 N 2 ^ 2 {^ + r) + p 2 XT-fa a a act kMN 2 ]T* n (6) 

+ / u 2 M(/ii + T)[r 2l a act {N 2 Xk - fi F j) - X 2 k 2 g t N 2 ] = 0, 
and the other variables are given by 



A% = 



J_ 

Xk 



a F n 2 + a a T f T- 
t^F + — 1 — —^g t T in 1 - — 



in 



9t ( 1 T in \ _ kg t ( T* n 



1 - -ir > v * = 1 



Oiact V M ) l a act V M 

1_ /ii+r ta l Mj' l2 -^ 1+T ) Im { 1 M 
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Depending on the particular values of parameters there can be between zero and two distinct steady 
states determined by the roots of equation It does not prove possible to determine stability 
of the steady state E^ in a closed form, and hence one has to compute the eigenvalues of the 
corresponding Jacobian numerically. 

In a similar way, when A\ ^ and A 2 = 0, there are again two options. The first one, denoted 
by E*i, describes the case when T* n = T* = T 2 * = 0, and the other variables can be found as 



Xpik Xp\k \ XpikNi J Xpi \ XpikNi 

The steady state E^ is stable, provided the following conditions are satisfied 

p 1 Xkr 1 N 1 > p\Xkr 2 Ni + r^fip, a ac tPirilXkNi > p\g t X 2 k 2 Ni + a ac tru 2 fi F , 

(7) 

~iHFT\ {hf + 7) + P\XkN l [(up + 7) 2 + PfI ~ PiAfcJVi] > 0. 
The steady state E^ can undergo a steady state bifurcation when 

piXkrxNx = p\Xkr 2 Ni + rufip, or a act Pin'yXkNi = p\g t X 2 k 2 Ni + a^tn"/ 2 fi F , (8) 
and it can also undergo a Hopf bifurcation, when the system parameters satisfy the conditions 

p 1 Xkr 1 N 1 > plXkr 2 Ni + ri7//p, aactPm^XkNx > p\g t \ 2 k 2 Ni + a act ri^ 2 fj, F , 

(9) 

7^71 (/ip + 7) +piXkN 1 [(n F + 7) 2 + pfI ~ PiXkNi] = 0. 

Similarly to the case of Eg, the Hopf bifurcation of the steady state E^ takes place inside the 
flow-invariant manifold A 2 = Ti n = T\ = T 2 = 0. Analogously to the analysis of stability of the 
steady state Eg, one can note that for stability of the steady state En, the linear growth rate r\ 
should be sufficiently greater than the rate r 2 . 

The second option for the case A\ ^ and A\ = is E^ 2 , for which T* is different from zero 
and satisfies the quadratic equation 



-1* \2 
" in) 



gtl^act[ni{oiFm + a a r) + p 1 a a TN 1 Xk}(T* 

-gt[ri7 2 a act M(a F ^2 + a a T) - pjX 2 ^^^^ + t) + p^aaaactkMN^T^ (10) 

+fi 2 M(m + r)[ri7a aci (piA^iAA; - ^7) - p 1 X 2 k 2 g t N 1 ] = 0, 
with the other variables being given by 

7 



A\ 



piXk 



«F^2 + «gT / T in 

H -, j ^gUin 1 - -TT 

/X2(Ml+ r J V M 



F * = 9t A _ TL\ y* = kgt A T* n 



OLact \ M J 'yaact \ M 

9t T * f, _ T* n \ _ rg t 

lii+r in \ Mj' 2 ^i + ry m y M 
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Finally, when both A\ 7^ and A\ 7^ 0, there are again two options. The first one corresponds 
to a steady state E\ 3 with T* n = T* = T 2 * = 0, and the other variables being given by 

F * _ rir 2 -y [{piNi + N 2 )Xk - ^7] _ r x r 2 [(j?xNi + N 2 )Xk - fi F j] 

X 2 k 2 {r 1 N 2 +p 2 r 2 N 1 ) ' ~ X 2 k{r 1 N 2 +p 2 r 2 N 1 ) 

A * _ N 1 [XkN 2 (r 1 - pir 2 ) + HFir 2 p\] ^ _ N 2 [XkpiNi(pir 2 - ri) + /x^n] 
1_ Xk(nN 2 + pjr^) ' 2 ~ Afc(nJV 2 +^r 2 JVi) 

Numerical computation of the eigenvalues of linearization shows that for sufficiently small values 
of 7, the steady state E[* 3 is stable; as 7 increases, stability is lost, but then regained again as 7 
increases further still. 

The final possibility corresponds to a steady state Ef 4 with the values of all variables being 
different from zero, so that T-* satisfies the quadratic equation 



gtictact [Xa a Tk(pir 2 Ni + p 2 r 1 N 2 ) +jrir 2 (a F fj, 2 + a a r)} (T* n ) 2 
[Xka a T(pir 2 Ni + p 2 riN 2 ) + jnr 2 (aFfJ-2 + ot a T)] 

(11) 

-X 2 k 2 fi 2 (fn + T)(plr 2 Ni + nN 2 )}T* n - n 2 M{n x + T)[X 2 k 2 gt(plr 2 Ni + n N 2 ) 

-Xk'ya.actrxriipxN-L + N 2 ) + n,Fl 2 a a ctnr 2 ] = 0, 

and the other variables can be found as 

Ni [n 2 M(m + T)(ri7a aci - p\Xkg t ) + ja a a ac tTg t T in (T* n - M)\ 



a; 



-4; 



M^27«ac*n(/Xl +t) 

N 2 [/z 2 Qm + T){r 2l a act M + Xkg t {T* n - M)) + p 2l a a a act Tg t T* n {T* n - M)] 

M fj, 2 ja a ctr 2 (jii + r) 

Si ^ y* = k 9t L T* n 



a-act V M J l a act V M 



9t - T *n T 2 = J 9 *, ; T*Jl TL 



fxx+r m \ M ; ' 2 M 2 (Mi + r) m V M 



To illustrate how stability of different steady states is affected by the system parameters, we 
show in Fig. [1] regions in the k-X parameter plane where different steady states are stable depending 
on whether infection and autoimmunity affect the same of different cell populations. This figure 
suggests that for sufficiently small rate of infection A, the only possible stable steady state is the 
disease- free state E%, in which cell populations A\, A 2 and T in are maintained at their respective 
carrying capacities, and there no free virus or infection, thus implying the absence of activated T 
cells T\ and T 2 . As the infection rate A increases, this disease- free steady state loses stability, and 
the system turns to a stable state of chronic infection E\^, where all cell populations are different 
from zero. When this steady states loses stability in the case of p\ = p 2 = 1 (i.e. when both cell 
populations A± and A 2 are the targets of both infection and autoimmunity), there appear to be no 
other stable steady state for any values of A and k, as illustrated in plot[TJa). At the same time, it 
is worth noting that some other stable steady states can be recovered when other parameters are 
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(a) (b) 




3456 3456 
k k 

Figure 1: Stability regions of different steady states depending on A and k. Parameter values are 
n = 0.2, r 2 = 0.1, Ni = 200, N 2 = 300, M = 300, g t = 0.2, a F = 0.0006, a act = 0.04, a a = 0.033, 
Hf = 1) A*i = 0.08, /X2 = 0.02, r = 0.001. (a) pi = P2 = 1- (b) pi = P2 = 0. Colour code denotes 
region of stability of a single steady state: Eg (violet), E*^ (green), E* 3 (red), E* (blue), and white 
denotes a parameter region where all steady states are unstable. 

changed. For instance, stability of the steady states Eg and En depends heavily on the relation 
between the two linear growth rates n and r 2 , so varying these may produce additional stable 
steady states. In the case pi = p 2 = 0, the range of possibilities is much wider, as shown in plot 
Hfb). It is possible to observe stable state E\q, where the cell population A\ is equal to zero with 
all other populations being positive. Another possibility is a state in which infection renders 
virus specific T cells ineffective. An important observation is that in all parameter regimes there 
appears to be a single stable steady state, hence no bi-stability is possible in the system. Whilst 
Figure Q] may not illustrate the complete catalogue of possible steady states, it provides an insight 
into how stability is affected by the system parameters. 

We summarize the above analysis of possible steady states of the system (pQ) in the following 
table. 
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Disease- free 
steady states 
V = F = 0, 
T in = Ti = T 2 = 


Ei 


all cell populations are equal to zero, 
always unstable 


E2 


Ti n = M, all other populations are 
equal to zero, always unstable 


E3 


Ai = N2, all other populations are 
equal to zero, always unstable 


Ei 


A2 = N2, Ti n = M, all other populations 
are equal to zero, always unstable 


E 5 


A\ = Ni, all other populations are 
equal to zero, always unstable 


Eq 


Ai = Ni, T^ = M, all other populations 
are equal to zero, always unstable 


Ef 


Ai = Ni, A2 = N2, all other populations 
are equal to zero, always unstable 


Es 


Ai = Ni, A 2 = N 2 , T in = M, all other 
populations are equal to zero, can be 
stable or unstable 


Steady states 
with a chronic 
viral infection 
V > 


Eg 


Ai = T^ = Ti = T 2 = 0, other 
populations are positive, can be 
stable or unstable 


E10 


Ai = 0, all other populations are 
positive, can be stable or unstable 


En 


A2 = T^ = T\ = T2 = 0, other 
populations are positive, can be 
stable or unstable 


E12 


A2 = 0, all other populations are 
positive, can be stable or unstable 


E13 


Tin = T\ = T2 = 0, other populations 
are positive, can be stable or unstable 


Eu 


all cell populations are positive, 
can be stable or unstable 



4 Numerical simulations 

In order to illustrate various dynamical regimes that can be exhibited our model, we solve the the 
system ([1]) numerically in different parameter regimes, taking account of results concerning stability 
and bifurcations of the steady states analysed in the previous section. Prior to performing simula- 
tions, we simplify the system by introducing the non-dimensional variables {A\, A 2 , F, T{ n , Ti,T2,Ai, V) 
and t as follows, 

t = g t t, A 1 = N X A 1 , A 2 = N 2 A 2 , F=(Ni + N 2 )F, 
T in = MTin, Ti=MT u T 2 = MT 2 , V = V V. 
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Substituting these variables into the system ([I]) yields 

dAi 
dt 

dA 2 



dt 

dF 
~dt 



r l A 1 (1 - Ax) - Pl \A x V - a a T 2 A u 
? 2 A 2 (l - 1 2 ) - \A 2 V - p 2 S a f 2 A 2 , 
PxNpM + (1 - N p )A 2 ] V-(jl F + apn + a a f 2 )F, 



= Tin ( 1 - T in ) ~ a ac tf in F, (12) 
dt \ ' 

-S- = a act f in F - f\(/Ii + f), 
dt 

dT 2 ^ ^ ~ 
= tTi - [i 2 T 2 , 

dt 



Z = kF-jV, 
dt 

where the modified parameters are given by 



^ ri ^ r 2 t XV ^ 7 t k(Nx + N 2 ) 
n = — , r 2 = —, A = , 7 = —, k 



9t 9t 9t 9t 9tV Q 

a a M ^ a F M ^ a act (iVi + N 2 ) ^ r 

OL a = , «F= , CXact = , T = — , 

9t 9t 9t 9t 

/Ui _ fi 2 hp 

= — > = — , = — , 

9t 9t 9t 

and we have introduced one additional parameter 



N, 



p Nx+N 2 

The above non-dimensionalization reduces the number of free parameters by three, thus reducing 
the overall complexity and simplifying the survey of parameter space. To simplify the notation, we 
will drop hats for variables and parameters in the system (|12p . The values of system parameters 
used in simulations of system (|12p are the rescaled values of the parameters used in Fig. [TJ 

4.1 Single infection 

First, we consider the situation when the host experiences a single infection by a virus, and one 
is interested in the subsequent dynamics of the immune response against this infection, as well as 
possible autoimmune reaction. Following Wodarz et al. [52j[53] and Vickers et al. [50], we define a 
certain threshold value, below which the infection is considered extinct; for simulations presented 
below this threshold was chosen to be 10 -8 and applied to the number of infected cells F. Initial con- 
dition was taken to be {A^O), A 2 (0), F(0),T in (0), Ti(0), T 2 (0), V(0)) = (0.9, 0.0333, 0, 0.9, 0, 0, 0.05) 
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Time (a.u.) Time (a.u.) 



Figure 2: Numerical solution of the system (|12p for p\ = P2 = 1- Parameter values are n = 1, 
r 2 = 0.5, iVp = 0.4, g t = 0.2, a a = 50, \i F = 5, a aci = 100, a F = 1, m = 0.4, /x 2 = 0.1, r = 0.005, 
k = 105, 7 = 10. (a) and (b) A = 0.1, (c) and (d) A = 0.5, (e) and (f) A = 1, (g) and (h) A = 3. In 
all plots, the colours denote infected cells F (blue), activated T cells T\ (red), and autoreactive T 
cells T 2 (black). Time is measured in arbitrary units (a.u.). 
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for all simulations, representing the fact that before the infection there are no infected cells and no 
activated T cells. 

Figure [2] shows the dynamics of system (|12|) for the case p\ = P2 = 1, which corresponds to 
a situation when both types of cells A\ and A2 are targets of both infection and autoimmunity. 
For sufficiently small values of infection rate A, the infection is being completely cleared: after the 
initial peak, the number of infected cells is monotonically decreasing, and the system approaches a 
stable steady state Eg, as shown in plots (a) and (b). This is the case of a normal disease clearance, 
where immune response of the host is able to completely clear the infection without causing an 
autoimmune reaction. Due to the exhaustion of the pool of infected cells, the activation of naive T 
cells stops, and the population of activated T cells is then slowly diminishing. The same happens to 
autoreactive T cells, whose number reaches its peak slightly later than the population of activated 
T cells, as these cells have to be derived from the population of regular activated T cells T\. It 
is worth noting that unlike an earlier model of Iwami et al. [23, 24J, the present model is able to 
support the initial viral expansion and subsequent clearance of infection by the immune system. 

For higher values of A, the disease- free state J5| loses stability, and the system evolves to a 
stable equilibrium J5* 4 , which describes the state of chronic infection. In this case, the immune 
system of the host is unable to clear the infection, and as a result it persists at a constant level. As 
A increases further, the steady state £7* 4 loses stability via Hopf bifurcation, giving rise to stable 
periodic oscillations, as shown in plots (e) and (f). In this case, one observes episodes of high 
viral production (relapses) with long periods of quiescence (remissions). Such dynamics have been 
observed in a number of autoimmune diseases, such as MS, autoimmune thyroid disease, uveitis etc. 
[Hll5tl39]. An important note here is that none of the subsequent reactivations of the virus requires 
any exogenous factors, but rather the system itself cycles through periods of relative quiescence 
and viral release. During each viral episode, the number of infected cells increases, and this triggers 
rapid activation of naive T cells, which in turn suppresses viral production, leading to a decrease 
in the number of infected cells. As the autoreactive T cells in the simulation shown in Fig. [2] have 
a much longer lifetime, their number decreases much more slowly during the periods of relative 
suppression of infection. For even higher values of A, this periodic solution becomes unstable, and 
instead the system evolves into another periodic orbit arising from the Hopf bifurcation of the 
steady state E* 2 , in which case the second population of target cells Ai goes to zero. 

Next, we consider another biologically plausible scenario when cells A2 are only targeted by 
infection, and cells A\ are only affected by autoimmunity, which is described by Pi = P2 = 0. 
Figure [3] illustrates the dynamics of the system (]12p in this case. For sufficiently small values of 
A, one observes normal clearance of infection similar to the case p\ = P2 = 1, and as A increases, 
the system tends to a stable state of constant chronic infection E*^, as shown in plots (a) and (b). 
When the steady state E\^ loses stability, rather than develop sustained oscillations as in the case 
Pi = Vi = 1) now the system goes instead to a stable steady state E* 3 shown in plots (c) and (d), 
which has Tj n = T\ = T2 = 0. This situation describes a state, in which the numbers of naive, 
activated and autoreactive T cells are all zero. While this is not biologically realistic, functionally 
it resembles exhaustion in which virus specific T cells are rendered ineffective and therefore the 
effective population size is reduced to zero. For higher values of A, there is another stability switch, 
and the system evolves toward a stable equilibrium E* . This behaviour is shown in plots (e) and 
(f), and it describes a situation when the first population of target cells A\ goes to zero. When 
A is increased further still, the steady state E* becomes unstable, and one observes a state of 
autoimmunity represented by stable periodic solution arising from a Hopf bifurcation of a chronic 
steady state E*^. 

In the case when only one of p\ and P2 is different from zero, the system can exhibit behaviours 
similar to the cases described above, with transitions between stable steady states and periodic 
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c 1 0" 3 (a) (b) 




Time (a.u.) Time (a.u.) 



Figure 3: Numerical solution of the system (I12p . Parameter values are the same as in Fig. [2j except 
for Pl = P2 = 0. (a) and (b) A = 0.9, (c) and (d) A = 1.5, (e) and (f) A = 5, (g) and (h) A = 9. In 
all plots, the colours denote infected cells F (blue), activated T cells T\ (red), and autoreactive T 
cells T2 (black). Time is measured in arbitrary units (a.u.). 
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Figure 4: Numerical solution of the system (|12p . Parameter values are the same as in Fig. [21 except 
for p\ = p2 = 0. (a) and (b) A = 0.2, (c) and (d) A = 1. In all plots, the colours denote infected 
cells F (blue), activated T cells T\ (red), and autoreactive T cells T2 (black). Time is measured in 
arbitrary units (a.u.). 



solutions of different origins. Numerical simulations suggest that in most cases, for sufficiently 
large values of the infection rate A, the system approaches a stable periodic orbit arising from the 
Hopf bifurcation of the steady state E* A . Such periodic solution corresponds to the above-mentioned 
state of autoimmunity with remissions and relapses. 



4.2 Multiple infections 

Since for many viral infections it is realistic to expect subsequent exposures of a person to the 
same virus, we now consider a scenario where someone who has recovered from a primary infection 
or currently has a chronic viral infection experiences a secondary viral challenge with the same 
virus. It is known that the timing of secondary infection plays an important role in determining 
the progress of infection, as well as the immune dynamics [341 147j . 

Figure |4] illustrates the dynamics during a secondary viral exposure in the parameter regime 
with normal clearance of infection and during the chronic infection. In the case of normal clearance 
shown in plots (a) and (b), one can observe that due to a much slower decay of activated T cells, they 
still remain at a non-negligible level following the primary infection. This means that the second 
infection produces a significantly smaller number of infected cells. The later a secondary viral 



15 



challenge occurs, the higher will be the resulting number of infected cells, and correspondingly the 
higher will be the numbers of activated and autoreactive T cells during a secondary infection. We 
find that the number of T 2 cells during a secondary infection is often higher than during a primary 
infection, such that the exacerbation of the autoimmune reaction is relatively much greater than 
the antiviral response, as shown in plot (b). The actual level of T\ and T 2 cells depend on the timing 
of secondary infection: if a secondary infection takes place sufficiently close in time to the original 
infection, the number of T\ cells will only increase by a very small amount, while the number of 
T 2 cells will exceed that number during the primary infection. If, however, a greater time elapses 
between the two infections, the number of activated T cells T\ will be slightly greater, but still it 
will never exceed the number of T\ cells during a primary episode. At the same time, the number 
of cross-reactive T cells T 2 in this case will be lower than during the primary infection. 

Plots (c) and (d) show that when the system is chronically infected, due to the significant 
amount of activated T cells, the secondary infection does not lead to a major increase in the 
number of infected cells, and as a result the infection is quickly cleared to the same chronic level 
as before the secondary infection. We have also analysed the influence of secondary exposures on 
the dynamics of recurrent infections, and in this case the effect of subsequent viral challenges is 
quite small in that it does not change the amplitude or period of oscillations except producing a 
small additional peak in the number of infected cells immediately after the infection. From this we 
conclude that in the parameter regime when the infection is recurrent, the main role is played by 
the periodic nature of the system itself, and it is this that causes relapses and remissions, rather 
than subsequent infections. 

4.3 Treatment 

As a next step in the analysis of onset and dynamics of autoimmunity, it is instructive to consider 
a practically important issue of therapeutic intervention. Numerical simulations presented earlier 
suggest that in the case when autoimmunity is triggered by infection with a virus, there is no 
further need for infection to maintain the periodic state of relapses and remissions. Let us now 
consider how the autoimmune dynamics changes upon the introduction of a therapy, aimed at 
reducing the number of autoreactive T cells. It will be assumed that such treatment can have two 
potential impacts: it can either act by eliminating the autoreactive T cells, thereby reducing the 
overall burden of autoimmunity, or it can transform those autoreactive T cells into cells with a 
higher activation threshold, i.e. activated T cells. These scenarios correspond loosely to treatments 
with therapeutic monoclonal antibodies that delete populations of cells [13j . or with drugs such as 
cyclosporin that inhibit the activation of T cells |31j . In clinical studies, these therapies cannot 
yet be made to selectively target autoreactive cells, which remains an important ongoing goal 
of treatment. Mathematically, we model such an intervention by modifying the fifth and sixth 
equations of the system (fT]) (and, correspondingly, (fT2l) ) as follows 

a act T in F - T1O1 + t) + p 3 CT 2 6(t - t ), 

(13) 

r^- ^ 2 T 2 -CT 2 6(t-t ), 

where C is the rate at which autoreactive T 2 cells are destroyed by treatment, p% shows the fraction 
of autoreactive T cells that can be converted into activated T cells (when p 3 = 0, the treatment 
is only responsible for reducing the number of T 2 cells), #(•) is the Heaviside function, and to is 
the time when the treatment is introduced. We assume that initially a person is exposed to a viral 
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Figure 5: Temporal dynamics during treatment of an autoimmune state. Parameter values are the 
same as in Fig. [21 and p$ = 0. (a) and (b) p\ = p2 = 0, A = 10, C = 30; (c) and (d) p% = 1, pi = 0, 
A = 0.65, C = 2; (e) and (f) p\ = 0, p2 = 1, A = 0.9, C = 2. Time is measured in arbitrary units 
(a.u.). 
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Figure 6: Temporal dynamics during treatment of a chronic infection. Parameter values are the 
same as in Fig. [2J pi = p2 = pj, = 0, A = 1, and C = 30. Time is measured in arbitrary units (a.u.). 



infection and develops some sort of autoimmune response, and in response to this a treatment is 
introduced in order to reduce or eliminate autoimmunity. 

Figured] shows the effects of treatment on the dynamics of immune response in the regime, where 
before treatment the system exhibited sustained periodic oscillations corresponding to relapses and 
remission during autoimmunity. All the simulations were performed for the case P3 = 0, but it 
has been checked that the same results are obtained for p% > 0. In the case when p% = p2 = 
(i.e. cells A<± are the only target of infection, and cells A\ are only affected by autoimmunity), once 
the treatment is introduced, the number of autoreactive T cells reduces significantly, as shown in 
plots (a) and (b). Although the oscillations persist after treatment, the number of A\ is greatly 
increased, and for sufficiently high rate of treatment C, it stays very close to 1. For p\ = 1 and 
P2 = (infection targets both A\ and A2 cells, but only A\ cells are affected by autoimmune 
response), plots (c) and (d) indicate that the introduction of treatment leads to a stable steady 
state -EJ^, in which all T cells (Tj n , T\ and T2) are eliminated. It is noteworthy that whilst the 
autoimmune reaction is eliminated in this scenario, this also leads to an increased level of persistent 
infection. When p\ = and P2 = 1 (only A2 cells are a target of infection, and both A% and A2 cells 
are affected by autoimmunity), treatment leads to suppression of oscillations and establishment of 
a stable chronic state E*^, as illustrated in plots (e) and (f). Due to a high rate a a , at which 
autoreactive T cells T2 destroy infected cells, once the population of these cells is reduced, the 
resulting chronic state E14 s characterized by a higher level of infected cells. In the case when 
Pi = P2 = 1) behaviour of the system under treatment is qualitatively similar to the case of p\ = 1 
and P2 = 0, when all T cells are eliminated. 

In Fig. [6] we illustrate how treatment affects the dynamics in the case of chronic infection. One 
can observe that treatment leads to a significant reduction in the number of autoreactive T cells, 
prompting a substantial increase in the number of A\ cells. However, treatment does not completely 
eliminate the infection which remains chronic, and similarly to the treatment of an autoimmune 
state, while the treatment reduces the level of autoimmune reaction, it simultaneously leads to a 
relative increase in the number of infected cells. This is biologically reasonable since inhibiting 
autoimmunity can lead to the reactivation of infection [21] . 
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5 Discussion 



In this paper we have developed and analysed a mathematical model for the dynamics of the immune 
response against a viral infection and the associated onset of autoimmunity. Having introduced 
separate populations of target cells that can be affected by infection and/or autoimmunity, as 
well as different compartments for T cells with different activation thresholds, we have studied 
how the outcome of the immune response can be the clearance of infection, chronic infection or 
recurrent infection. In the case of normal clearance, there are no lasting immune consequences 
for the organism. The number of autoimmune T cells is small and decreases with time after the 
clearance of infection. Chronic infections are characterized by a constant level of activated and 
autoreactive immune cells, which keep infection in check but do not clear it. Another possible 
scenario for persistent infection is when virus specific T cells is exhausted, and hence the effective 
T cell population is reduced to zero. Recurrent infections are characterized by inability of the 
immune response to clear an infection, which results in relatively long periods of remission followed 
by relapses with a substantial viral production and a large number of cells being infected. In this 
study we have focussed on chronic infection driven by an increase in the infection rate A. Studies in 
other models show that chronicity can also be influenced by other factors such as the frequency of 
antigen specific T cells [35] . During recurrent infection we observe higher numbers of autoimmune 
cells (compared to the case of chronic infection), which exacerbate the relapses and cause more 
damage. Notably, there is a difference in the outcome of the infection, depending on whether 
infection and autoimmunity affect the same or different organs. 

Numerical solutions reveal a number of potentially important features in the dynamics of the 
autoreactive cells. When a person experiences a secondary challenge with the same virus, the 
model shows that in the case of normal clearance the number of cells getting infected is smaller 
in subsequent infections due to a limited number of T cells remaining activated at the time of 
secondary infection. Although in this case the infection is still successfully cleared, in a number 
of circumstances the second and subsequent peaks in autoreactive cell numbers exceeded that seen 
during the initial response to virus (Fig. 2] (b)), thus causing an additional immune challenge. 
This could implicate repeated cycles of T cell expansion, or multiple linked rounds of infection as 
crucially underpinning the development of frank autoimmune disease. This would be consistent 
with multi-hit models of autoimmunity (e.g. [3\) and also with circumstances in which more than 
one infectious episode is necessary to precipitate frank disease. In contrast, secondary infection 
during a chronic infection does not have a significant effect on the dynamics. The same conclusion 
holds for the recurrent infections, thus indicating that in this case it is the intrinsic dynamics of the 
interactions between the immune system and the virus that causes remissions and relapses rather 
than the fact that a person experiences further infections. 

In the form of the model where the population of infected cells is different from that which is the 
target of autoimmunity, we found that under some conditions, low levels of viral persistence could 
be associated with high levels of autoreactivity (Fig. [3] (e)&(f)). We also demonstrated that the 
system can approach a state that resembles T cell exhaustion (Fig. [3] (c)&(d)). This is intriguing 
because inactivation of the immune response in the face of chronic viral infection and in tumours 
is a well described and important area of ongoing investigation [26] . 

We have also studied the effects of treatment aimed at reducing the number of autoreactive 
cells on the dynamics of the immune response. In the case of the recurrent autoimmune state, 
when infection and autoimmunity occur in different organs, such treatment leads to a substantial 
improvement of the situation, significantly reducing the number of autoreactive T cells. It is impor- 
tant to note that although this does not eliminate episodes of relapses /remissions, they have a much 
less prominent impact on the numbers of susceptible cells and hence cause significantly less damage 
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than the full-blown autoimmunity. When only one of the susceptible cell populations is the target 
of infection but both of them can be affected by autoimmunity, introduction of treatment leads 
to suppression of relapse/remission oscillations and establishment of a state of chronic infection. 
When initially the infection is chronic, treatment does not qualitatively change the dynamic state 
of the system, but leads to a reduction in the number of autoreactive T cells. One should note, 
however, that due to reliance of immune system on autoreactive T cells to contribute to control of 
infection, any treatment aimed solely at reducing the number of these cells can inadvertently lead 
to a higher level of persistent chronic infection. This trade-off between the ability to control the 
infection and at the same time to minimize the undesired effects of autoimmunity is somewhat sim- 
ilar to the problem in chemotherapy where an effective treatment of a tumor may have a negative 
impact on the overall immune status, thus requiring some sort of adjuvant therapy. At the same 
time, if the role played by autoreactive T cells in the clearance of infection is not so significant, the 
above problem becomes less serious, and it is possible to achieve efficient control of autoimmunity 
without compromising host's ability to fight infection. 

There are a number of interesting potential extensions to this work which may be possible in the 
future. While we started with the concept of tunable activation thresholds, the implementation we 
chose here was deliberately simplified to ensure the tractability of the model. Previous theoretical 
and experimental work [H [21 [181 021 EH SS] has stressed the dynamic nature of the tuning process 
which can be observed over timescales of minutes to days [431 [48J . It would therefore be interesting 
to explore the dynamics of this process and the effects that different parameters for the kinetics of 
tuning might have on the development of autoreactive T cell populations. Another feature of the 
current model is that the autoimmune response does not persist when the virus is cleared. Since 
there is good experimental evidence that autoimmune responses can be self-sustaining and chronic 
[25], developing the model to explore this behaviour will be an important further development of 
the model. 
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